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Abstract. We show how to adapt the ideas of local energy and momentum conservation 
in order to derive modifications to the Gross-Pitaevskii equation which can be used 
phenomenologically to describe irreversible effects in a Bose-Einstein condensate. Our 
approach involves the derivation of a simplified quantum kinetic theory, in which all processes 
are treated locally. It is shown that this kinetic theory can then be transformed into a number of 
phase-space representations, of which the Wigner function description, although approximate, 
is shown to be the most advantageous. In this description, the quantum kinetic master equation 
takes the form of a Gross-Pitaevskii equation with noise and damping added according to 
a well-defined prescription — an equation we call the stochastic Gross-Pitaevskii equation. 
From this, a very simplified description we call the phenomenological growth equation can 
be derived. We use this equation to study i) the nucleation and growth of vortex lattices, and ii) 
nonlinear losses in a hydrogen condensate, which it is shown can lead to a curious instability 
phenomenon. 
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1. Introduction 

There is a growing consensus that in some sense the dynamics of a trapped Bose-Einstein 
condensate can be properly described by the Gross-Pitaevskii equation, provided appropriate 
statistical assumptions are introduced. This was first suggested by Svistunov [[j], ^, gp, an d 
numerical experiments were carried out by Damle and Sachdev Marshall, New and 
Burnett [||, and most recently by Davis et al All of these considered the Gross-Pitaevskii 
equation for a Bose gas inside a box, with random initial conditions, such that the total energy 
and number of particles would correspond to those of a partially condensed system in thermal 
equilibrium. The numerical experiments then verified that this initial ensemble developed with 
time into a thermodynamic equilibrium ensemble as a result of the time evolution induced on 
each member of the ensemble by the Gross-Pitaevskii equation. The recent work of Davis 
et al ^ also checked the dynamics of the resulting equilibrium state against Morgan's [Q] 
theoretical predictions from quantum field theory, and found that these were in agreement 
with each other. 

The work of Steel et al. took a different, and more rigorous view, that the correct 
physical relationship between the condensate wavefunction with random initial conditions and 
the physical field operator is made by means of a truncated Wigner function representation. It 
is known ^ that for the standard Bose-Einstein condensate Hamiltonian, the Wigner function 
obeys a generalized Fokker-Planck equation, in which derivatives up to third order occur. 
However, if the amplitudes of the Wigner function phase-space variables are large, the third 
order derivative terms are small, and can be neglected — this "truncation" gives the method 
its name. However, we also find that the expected second-order derivative terms, which 
correspond to noise, are for this Hamiltonian exactly zero, and thus what remains corresponds 
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to a Liouville equation for an ensemble of wavefunctions which obey the Gross-Pitaevskii 
equation. 

However the description can only be valid provided the amplitudes can be regarded 
as large, and this must be the case for all modes. Clearly, the higher energy (or smaller 
wavelength) modes will always have a very small amplitude, and cannot be described 
by a truncated Wigner function. Indeed, if no truncation is done, then the result of 
the quantum mechanical constraints — essentially Heisenberg's uncertainty principle — on the 
Wigner function leads to a delta function spatial correlation function, corresponding to 
infinite fluctuations at each point in space, which cannot be simulated. In all simulations, 
however, the implementation of the spatial differentiation requires a spatial grid, whose 
spacing corresponds to the inverse of the highest spatial frequency used, thus providing a 
cutoff in both space and energy. The results must then be cutoff dependent. 

The work of Davis et al ^ is the first to put down explicitly that a cutoff in energy 
or spatial frequency must be imposed on this Gross-Pitaevskii equation equation. However, 
the basic idea, that high and low energies require different treatments is that which forms the 
basis of our quantum kinetic theory, which gives a formalism for coupling together a fully 
thermalized noncondensate band, consisting of particle with energies more than a certain 
value Er, to a condensate band, consisting of all lower energy particles. The non-condensate 
band can be treated by the quantum Boltzmann equation, while the condensate band should 
in principle be treated fully quantum mechanically. 

In this paper we will combine the Wigner function ideas of Steel et al. [[| with quantum 
kinetic theory, and as well introduce the idea of local energy conservation as suggested by 
Zaremba et al. [10 TTJ] to derive a method of coupling the thermalized noncondensate band to 
the Gross-Pitaevskii equation. This results in a Gross-Pitaevskii equation with added noise — a 
stochastic Gross-Pitaevskii equation. 

The paper consists of three main parts: Sect.[| implements the idea of local energy and 
momentum conservation in quantum kinetic theory and results in a relatively simple master 
equation for the interaction of a condensate with a fixed bath of non-condensed atoms; Sect.|] 
shows how to implement the Wigner function method to transform the master equation into 
a c-number equation with random initial conditions, and to compare with P- and Q-function 
methods; Sect develops a very much simplified phenomenological equation, which we show 
can be applied to the two-body loss problem in the hy drog en condensate experiments, and 



to the stabilization of quantized vortex arrays. In Sect. 6.1 we adapt the phenomenological 



growth equation to include losses from dipolar relaxation in a hydrogen condensate, and show 
that the nonlinear losses so introduced could lead to a "boom and bust" instability in which 
the condensate grows and collapses repeatedly, although the conditions under which hydrogen 
condensates are presently formed make it difficult to say whether this could be observed in 
practice. 



2. Application of local energy conservation to quantum kinetic theory 

In a formulation of the kinetic theory of non-condensed vapour in interaction with a 
condensate, Zaremba, Nikuni and Griffin [[ll]] have used the concepts of local energy and 
momentum conservation to develop appropriate equations of motion for the system in a 
formulation based on Hartree-Fock-Popov methodology [ |l2[ |l3|, |l4[ [l5[ ]. This is in contrast 
to our own 

jT|, HI HI IB H' H formulation of quantum kinetic theory, in which 
energy conservation is expressed in terms of transitions between eigenstates of the condensate. 
While it is clear that a description in terms of eigenfunctions must be more accurate, it yields 
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equations which are not easy to handle exactly, and thus the greater accuracy can in practice 
be an illusion. 

The full formulation of our quantum kinetic theory is to be found in QKV ||| which 
depends strongly on QKIII. The major improvement in QKV is a full consideration of the 
mean field effects of condensate on vapour and conversely, but the general methodology is 
very similar. In these we divide the excitation spectrum at an energy Er, above which the 
excitations can be treated as particle like. The field operator is written (in the Schrodinger 
picture) as 

i/)(x) =^ NC (x)+0(x) (1) 

where the first operator represents the modes with excitation energies above Er (and is called 
the noncondensate band field operator), while 0(x) represents the excitations with energies 
below Er (and is called the condensate band field operator). This division is chosen and fixed 
for any given problem, and is to some extent arbitrary, since it is expected that the higher 
modes of the condensate band will be very close to particle-like, and in practice will also be 
thermalized. 

Let us look at the derivation of energy conservation used in QKIII, using the notation 
used there. On p539 we examine a term in the master equation involving Hj C as defined 
in (QKIII. 8), which contains one 0(x) or 0^(x), where 0(x) is the field operator for the 
condensate band. Explicitly, the term in Hamiltonian can be written 

H ?C = J d 3 xZ 3 (x)^(x)+h.c. (2) 
In this equation we have defined a notation 

Z 3 (x) = uVn C (x)V'nc(x)'i/>nc(x). (3) 
Substituting into the master equation (QKIII. 20), terms arise which are of the form 

--L Jd 3 xJ d 3 x' J dTTr{z 3 (x)Zl(x.',-T)p NC } 

x^{x)d,(yf,-T)pc{t). (4) 
Here we have used the notation 

Z 3 (x,t) = e * H *ct/h z ^ x y- t H NC t/h (5) 

0(x,*) = e tHot/h cj)(x)e- tHot/R . (6) 

In QKIII we expanded the operators </>(x) in eigenoperators of the condensate band 
Hamiltonian Hq, so as to be able to perform the integral over the time t, and thus arrive 
at the final master equation, which would not involve <fr operators at different times. This 
method does not require that there be a condensate. If we assume that there is a condensate 
present, we use this knowledge to make an estimate of the time development in the time r. 



2.1. Local energy conservation 

The basis for this concept is given by the works of Zaremba, Nikuni, Griffin, and co-workers 
], which themselves draw on the work of Kirkpatrick and Dorfman [[27j ^8], p9| |, 
and can be seen phenomenologically by use of the density-phase description of the Gross- 
Pitaevskii equation 

ih ^§T^ = -|^V 2 e(x,t) + y T (x)C(x,t)+ W |^(x,t)|^(x,i). (7) 

| In this p ape r will use the notation: QKI fo7^, QKII for QKIII for QKIV for l^, QKV for 
QKVI for and QKVII for ^ 
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By writing now 



£(x, t) = y/n c (x, t) exp[ze(x, t)} (8) 
the Gross-Pitaevskii equation takes the form 

dnc(x, t) 



8t 



= -V - [v c (x,i)n c (x,i)], (9) 



<99(x,i) , , 1 . 

h m = -Mc(x i t)--m t ; c (x,t) 2 (10) 

= -e c (x,t). (11) 

in which the condensate velocity vc(x, t), and the local chemical potential ^tc( x , t) are given 
by 

v c (x,t)E-Ve(x,(), (12) 
m 



A*c(x,t) = v °\' ; +y T (x)+^n c (x,t). (13) 

This transcription of the Gross-Pitaevskii equation is the appropriate basis for a description 
of condensate behaviour when the condensate density tiq (x, t) is slowly varying in space and 
time, so that we may consider the dependence on space and time to arise largely from the 
phase 9(x, t). This leads to the hydrodynamic approximation in which the Laplacian term in 
dTl) is dropped and thence to the energy conservation equation 



dt 



Vt( x ) + imt;c( x , t) 2 +un c (x, t) 



(14) 



The static solution of the hydrodynamic description is the Thomas-Fermi description of the 
condensate wavefunction, in which 

vc(x,i)=0, (15) 
Mc(x,i)=M, (16) 

, c (x,*) = ^M, (17) 

u 

e(x,t) = (is) 

In the hydrodynamic regime, the local chemical potential /ic (x, t) is a strictly local quantity, 
and is a simple function of the density rtc(x, t). The local energy density ec(x, t), defined 
in ([TT]) is explicitly equal to the derivative of the energy density (also evaluated using the 
hydrodynamic approximation) 

E(x,t) = imw c (x,t) 2 n c (x,i) + ^ T (x)7i c (x,t) + in c (x,<) 2 (19) 

with respect to nc(x, t). Thus, when considering the possible addition of a particle to the 
condensate, it is intuitively appealing to say that this takes place at a position x, and that the 
energy added to the condensate by this process is e(x, t). 

This is given some further backing by the knowledge that energy conservation in 
quantum mechanics depends on frequency matching, and that to the extent that nc(x, t) and 
e(x, t) depend slowly on time, the relevant frequency will be e(x, t). If nc(x, t) and vc(x, t) 
are both slowly dependent on space, the wavelength of the wavefunction will correspond to a 
momentum tovc (x, t), and thus momentum conservation will involve this momentum. 
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2.2. Application to quantum kinetic theory 

The major task in developing a quantum stochastic description is to find a simple way of 
expressing the term <^>(x', — r) in terms of <fi(x, 0), and we would like to follow Zaremba et al 
[ |ll| ] in doing this. We assume that the evolution over this short time can be approximated 
by simply the phase evolution of the condensate according to ( |Tol|lTl ). We also want to 
express the result in terms of operators at the same location, which we take as the midpoint 
u = (x + x')/2. Thus, using the symbol y = x — x', we write 

^(x)^(x',-r) «0t( u )^ u ) eX p{i[-e(y/2,t) + e(-y/2,t-r)]} (20) 
~ f (u)(/)(u)exp j~[-mv c (iM) • y - e c (u,t)r]| . (21) 

There are two main assumption here: 

i) There is an assumption that the state of the condensate band is dominated by a single 
condensate wavefunction. However, the wavefunction may have any form, and of 
course may also be time dependent. It is possible that a more accurate method could 
be developed using hydrodynamic quantization, such as that developed by Marques 
and Bagnato [|30|], which yields operator equations almost the same as the standard 
hydrodynamic formalism outlined above. 

ii) The approximation (|TJ) requires that only small y and r contribute. This will 

be valid when the function Tr lZz(x)Z^(x!, — t)pnc|> which is convoluted with 

(/)^(x)(/>(x', — t) in (Q), is almost local in space and r. This kind of behaviour is expected 
if the noncondensate band is fully thermalized at a sufficiently high temperature. 

Both of these assumptions are used by Zaremba et al JTT[ | 

2.3. Final form of the master equation 

By carrying out the procedures used in QKIII, we can finally get the master equation: 
Pc(t) = J JdPx f (x) - y T (x) - 2 Mj0NC (x,i) - i^t(x)0(x)^ 0(x) , Pc 

(22) 

+ /d 3 x(GW[x,e c (x,t)] (20(x) / j c t (x)- / 9G0 t (x)^(x)-0t(x)0(x)pc) (23) 



+ G(->[x,ec(x,t)] (20t(x)p c ^(x)-pc^(x)0 t (x)-0(x)^( x ) pc ) (24) 

+ M[x] (2C/(x)p c C/t( x ) - p c £/t(x)[/(x) - U\x)U{x)p c ) (25 ) 

+ £7< + >[x,ec(x,t)] (2V(x)p c ^ t (x)- / o c V t (x)y(x)-yt(x)y(x)pc) (26) 

+ £7<->[x, e c (x,i)] (2V+(x)pcV(x) - pcV(x)W(x) - V(x)V*(x)pc) ^ . (27) 

For compactness, we have defined 

C/(x) = t (x)0(x), (28) 

y(x) =0t( x ) ( / ) (x)^(x), (29) 

PNc(x,t)= /" rf 3 KF(K,x). (30) 
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Here F(K, x) is the phase space density of the noncondensate band in terms of wavenumber 
K and position x. 

The transition rates M, are defined by 

G (+) [x,e] = J d 3 Kx Jd 3 K 2 y , d 3 K 3 F(K 1 ,x)F(K 2 ,x)[F(K 3 ,x) + l] 

x^AK 123 -^f^)^(A Wl23 -|) (31) 

G^fx.e] = J* y d 3 Kx y"^K 2 y'd 3 K 3 [F(K 1 ,x) + l][F(K 2 ,x) + l] J F 1 (K 3 ,x) 

x^(AK 123 -^f^)^(A. 123 -|) (32) 

M(x) = -^2 y^Kx y d 3 K 2 F(K 1 ,x)[^(K 2 ,x) + l],5(K 1 -K 2 ) ( 5(cc )l -^ 2 ) (33) 

S(+)(x,e) = 2^ y d 3 K^(K llX )^ ( Wl - J) 5 (k, ( 3 4) 

£<->(x,e) = 2^ y d 3 Ki[F(Ki,x) + l]S ( Wl - |) * (k x - mVc fr'^ (35) 
In the above equations we use the notations 

hu(K,x)= T ^- + V T (x) (36) 

uji =cj(Ki,x) (37) 

AKi 23 ee Ki + K 2 K 3 (38) 

Awi 23 = uj\ + ll>2 - u>3. (39) 



2.3.1. Effect of the cutoff at Er The fact that 4>(x) has an upper energy cutoff means that 
the commutator [4>{x), ^(x')] is only an approximate delta function. This has no influence 
on the form of the master equation, but equations of motion derived from it for averages of <j> 
will automatically ensure that solutions remain in the correct subspace. 

2.3.2. Conservation of energy and momentum The choice made in Eqs (^,|T]) is very 
drastic, assigning as it does a single phase to all condensate band modes, and this has the 
consequence that the terms (^5]-|27|) are in fact zero. For the last two this is obvious, since the 
energy and momentum conservation delta functions in the definitions of E^' would require 
the equality of the energy and momentum of a particle in the condensate band with those of 
a particle in the noncondensate band, which is by definition not possible. The vanishing of 
the term (^5]) is not so immediately obvious, but momentum and energy conservation here 
mean that in fact this is a forward scattering term, which the methodology of QKIII and QKV 
explicitly removes from the irreversible part. 

The origin of this is the requirement that, at a given position x, there is a unique energy 
and momentum which a particle in the condensate band can have. This is actually only 
true for particles contained in the condensate itself, so that we can expect that there is a 
spread of values of ec(x, t) and Vc(x, t). As well as this, it should be borne in mind that 
the hydrodynamic approximation is the foundation of these results, and this is indeed an 
approximation. Inclusion of corrections to this approximation would also yield a spread in 
local momenta and energies available for particles in the condensate band. 



The stochastic Gross-Pitaevskii equation. 



7 



The effect of including a spread in values of ec( x , t) and vq(x, t) would be principally 
to give a nonzero value to the terms (fB]-|27|) in the master equation. The last two terms 
would certainly be very small, but the scattering term ( [25] ) should become significant. The 
growth terms ([22,23 ) would be modified, but not greatly. However the changes are only in the 
nonoperator coefficients — the basic structure of the master equation remains the same. The 
major point of simplification compared with the treatment of QKIII and QKV is the expression 
of the irreversible terms directly in terms of the field operators </>(x), which enables a much 
simpler analysis to be carried out. 



2.3.3. Time development equations for the noncondensate band The noncondensate band is 
described by the phase space density F(K, x), whose equation of motion is given in QKV, in 
the form of a quantum Boltzmann equations and with added terms to take account of transfer 
between the condensate and noncondensate bands. For simplicity, in this paper we assume 
that the noncondensate band is fully thermalized with temperature T and chemical potential 



3. Equivalent Wigner function treatment 

3.1. Stochastic differential equations for the Wigner function 

A master equation of the form we have produced can be treated using phase space 
representation methods [^]. These will give rise to equivalent stochastic differential equations 
for a c-number phase space variable a(x, t) whose averages are related to those of the operator 
0(x) using symmetric ordering: thus 

(a(x,t)> =(0(x)>, (40) 

(a*(x,i)) =(0 t (x)>, (41) 

<a*(x, iH x',t)) = ^(x)^(x)+^)^t( x)) ^ 

Using the operator correspondences in ^ Eq. (4.5.12), we can deduce the stochastic 
differential equation 

V 2 a(x, t) — Vr(x)a(x, t) — 2u/?nc(x, i)a(x, t) — it|a(x)| 2 a(x) > dt 

2m J 

+ {G (+ ) [x, e c (x, *)] - G<-> [x, e c (x, t)] - M[x] } a(x, t) dt 

+dW G (x,t) +ia(x)dW M (x,t). (43) 

Here the last two terms are Gaussian white noise terms. The quantity dW m ( x , t) is a real 
Wiener noise term, to be interpreted in the Ito sense [|3"T|], whose correlation functions are 

(dW M (x,t)}=0, (44) 

{dW M (x,t)dW M (x',t)) = 2M[x]5(x - x')dt, (45) 

while the other noise is complex, with correlation functions 

(dW G (x,t)} = (dW£{x,t))=0, (46) 

(dW G (x ) t)dW G (x',t))=0, (47) 

(dW G {^t)dW G ( X \t)) =0, (48) 

(dW G (x,t)dW G (x',t)) = S{X 2 X ° (GW[x,ec(x,i)] +G(-)[x,e G (x,t)]) dt. (49) 
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3. 1.1. Neglect of third order noise terms The equation is also approximate in that third order 
noise terms do arise along with the trilinear term. These terms arise as third order partial 
derivatives in the Fokker-Planck equation generated by the Wigner function representation of 
the master equation (p2|-[27|), and have no straightforward stochastic interpretation. When the 
variable a(x) is large, corresponding to large phase space occupation, these terms become 
small, as pointed out in [|J . 

3.1.2. Stoof's equation Stoof [[$2]] has developed a similar form of Wigner function 
stochastic Gross-Pitaevskii equation using path integral methods, and this has been used to 
analyse damping of condensate modes and the reversible formation of a condensate [[33], g4fl . 
The actual form of the damping is not the same as ours, but it is easy to see that to lowest 
order in the damping it is equivalent, and this is all we would claim for our derivation. The 
noise terms are equivalent to ours. 

In [ [33| ] a case is considered where one can take a "classical" limit of this Wigner 
approach — a low energy approximation and essentially always valid for the description of 
the dynamics of the condensate and its low-lying excitations — and in this approximation the 
P- Q- and Wigner function forms of the noise become the same. This leads to an approximate 
noise term which is of the the same form as that of the Q-function. 



3.2. Stochastic differential equations for P- and Q-f unctions 

3.2.1. Use of the P-function It is more conventional in quantum optics to use either the 
Glauber P-function or the positive P-function, instead of the Wigner function. The main 
advantage is that the approximation, necessary for the Wigner function, that third order noise 
terms must be neglected is not required — the resulting stochastic differential equations are 
exact. 

For either P-function form, we would replace the symmetrized product averaging rule 
with one involving normal products, so that we would obtain 

(a(x,t)> =Wx)>, (50) 

(a*(x,i)> =(^(x)), (51) 

(a*(x,t)a(x',t)) =<0t(x)0(x)>. (52) 

Using the operator correspondences in Eq. (4.5.9), we can deduce the same stochastic 

differential equation (^3j), but with the correlation functions 

(dW M (x,i)>=0, (53) 

{dW M (^,t)dW M {x',t)) = 2M[x]<5(x - x') dt, (54) 

(The same as for the Wigner function), while the other noise has the altered correlation 

functions 

(dW G (x,t)) = (dW G (x,t))=0, (55) 

(dW G (x,t)dW G (x',t)) = f «*(x) 2 , (56) 

111 

(dW G (x,t)dW G (x',t)) = --«*(x) 2 , (57) 

n 

(dW G (x,t)dW G (x',t)) = ,5(x-x')G (_) [x,e c (x,t)]dt. (58) 
The interpretation of (^3j) as a genuine stochastic differential equation requires that the matrix 
of noise coefficients 

G(->[x, ec (x,t)] -fa*(x) 2 \ 

fa(x) 2 G(-)[x,e c (x,i)] J PVj 
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should have only nonnegative eigenvalues. For higher temperature situations in which there 
is a substantial thermal component, this will certainly be true for all values of the variable 
a(x) which would turn up in a stochastic simulation. When this is not so, a Positive P- 
representation would be necessary. The experience of Drummond and co-workers [|5j ^] 
has shown this is in principle feasible, but application to experimentally realistic problems 
would be very difficult. 



3.2.2. Use of the Q-function In the case of the Q-function one should replace the 
symmetrized product averaging rule with one involving antinormal products, so that we would 
obtain 

(a(x,t)) = <<K X )>, (60) 

(a*(x,t)) =(^(x)), (61) 

(a*(x,t)a(x',t)) = (0(x)0t( x ')). (62 ) 

Using the operator correspondences in ^ Eq. (4.5.10), we can deduce again the same 
stochastic differential equation (j43|), but with the correlation functions 

(dW M (*,t))=0, (63) 
(dW M (x,t)dW M (x',t)) = 2Af[x]<5(x - x') dt, (64) 

(Again the same as for the Wigner function), while the other noise has the altered correlation 
functions 

(dW G (x,t)) = (dW G (x,t)) =0, (65) 
(dW G (x,t)dW G (x',t)) = -^*(x) 2 , (66) 

(dW G (x,t)dW G (x',t)) = ^a*(x) 2 , (67) 

n 

(dW G (x,t)dW G {x.',t)) = <5(x-x')G (+) [x,e c (x,i)]^. (68) 

The condition for this to be a genuine stochastic differential equation is that the matrix 

G(+)[x,e c (x,i)] fa*(x) 2 
-fa(x) 2 G(+)[x,e c (x,i)] 

Thus, although the Q-function always exists and is positive, this does not necessarily mean 
the stochastic differential equation is valid — see [ [3l| ] Ch.6. 



(69) 



3.3. Approximations and simplifications 

3.3.1. Effect of the cutoff These equations are expressed in a simplified form, since by 
definition the condensate band contains a finite range of energies, and hence of wavelengths. 

This means that the Laplacian, the delta functions, and the form of the trilinear term are 
modified by a projection into this band, as explained in QKIII and QKV. As well, the delta 
functions become delocalized, and the noise correlation functions are then well defined at 
x = x'. The actual implementation of this procedure is somewhat technical, but is essentially 
straightforward. 

A reasonable way to do this is to express the S(x — y) terms as a summation over a set 
of orthogonal condensate eigenfunctions. This amounts to truncating an exact expression 
for a delta function at the cutoff energy Er. But for sufficiently high energies the trap 
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eigenfunctions approach the harmonic oscillator eigenfunctions, so we may instead use the 
expression in terms of harmonic oscillator eigenfunctions y* (x) 

5(x-y)^5 fl (x,y) (70) 

= J2 £( x )w-(y) ( 71 > 

B(n)<B„ 

This representation necessarily gives an added noise which is only significant for x, y inside 
the region classically allowed for energy equal to Ep, and thus permits the use of fast Fourier 
transform methods. 



3.3.2. Effect of the cutoff on means and variances In application of the Wigner function to 
a field theory, the fact that each mode contributes an extra half of a quantum to averages — as 
shown in ( p2| ) — means that there will be a contribution proportional to the number of modes 
retained below the cutoff. For the Q-function the effect is similar, but the contribution is a 
whole quantum. 

Suppose M modes contribute to the summation in ([71]). If we define the particle number 
operators for the field theory and the various phase space representations 



N = I d 3 x(/) t (x)0(x) (72) 

N P = J d 3 xa P (x.)a P {x) (73) 

N w = J d 3 xa^(x)aw(x) (74) 

Nq = J d 3 xa^(x)a Q (x) (75) 

then the means and variances are related by 

(N P ) = (N) (76) 

(N w ) = (N) + M/2 (77) 

(N P ) = (N) + M (78) 

vax[N P ] = vax[N] - (N) (79) 

vax[N w ] = vax[N] + M/4 (80) 

vax[N Q ] = vax[N] + (N) + M (81) 



These results directly present the dilemma one is faced with in making a choice of 
representation. If we wish to consider only positive phase-space distributions, which therefore 
have a probabilistic interpretation, we can see immediately from these results that a P-function 
is only possible if there is at least a Poissonian number distribution — that is, the P-function 
does not always exist. (However the Positive P-representation can evade this restriction, at the 
cost of doubling the number of independent variables.) In contrast, there is no restriction on 
the Wigner and Q-functions, which are known to exist for all physical states. 

On the other hand, only the P-function has no "vacuum noise" contribution, proportional 
to M, the number of modes, unlike the Wigner and Q-functions. For these, the existence of the 
"vacuum noise" contribution means that a simple interpretation of the stochastic differential 
equation (^3|) as the equation of motion for the condensate wavefunction must require that the 
actual mean number of real particles (N) not be swamped by the vacuum contribution; that is 
we must require (N) ^> M at the very least. This means that the mean occupation per mode 
is very much greater than 1, and that a Bose-Einstein condensate must already be present. 
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Figure 1. a) The Wigner function for an N = 80 number state, and b) the corresponding 
cumulative Wigner function. 

For the P-function this is not a problem — M can be as large as we wish — even infinite if we 
desire. 

Finally, one must emphasize that the stochastic differential equation ( p3| ) is approximate 
for the Wigner function interpretation, even with the appropriate noise properties, because 
of the neglect of third-order noise terms, as previously noted. Fortunately, the condition 
that we van neglect these terms is the same as the condition that the vacuum modes not be 
dominant — that the occupation per mode is high for the modes of interest. On the other hand 
provided we have an initial distribution which is broader than the Poisson, and the temperature 
is sufficiently high for for the positivity condition of the noise matrix ( |59[ ) to be satisfied, the 
P-function interpretation of the stochastic differential equations ( pi3| ) is valid and exact. 



3.4. The low temperature case 

For sufficiently low temperatures it is clear that all the transition rates (|3"i|-|35|) are negligible, 
but the noise matrices ( |5^ , |69| ) no longer have positive eigenvalues. This leaves us with only 
two feasible choices — either a Positive P-function interpretation, which is numerically very 
difficult, or a Wigner function interpretation, which is approximate, but has the attractive 
property that the noise contributions all vanish, as can be seen from (@-0). This choice has 
been extensively investigated by Sinatra et al. J37|]. 



3.4.1. The stochastic Gross-Pitaevskii equation for a pure condensate To illustrate the 



procedure, we omit the noise and damping terms in (43), and consider the resulting time 



dependent Gross-Pitaevskii for a random function a(x, t). This takes the form 
i) a(x, t) satisfies the Gross-Pitaevskii equation: 

da h 2 

ih— = V 2 q + V T U)a + u\a\ 2 a. (82) 

at 2m 

The prescription of statistics in terms of symmetrically ordered products as in (0-0) 
means that the initial conditions must be random even in the case of an initial pure state. 
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There are two possibilities, associated with the fact that the total number of particles is 
now exactly conserved, because of the vanishing of the terms (|3l|-[3"5j). A straightforward 
interpretation would be to assume that the initial wavefunction of the condensate mode 
is £(x, 0) and use random initial conditions corresponding to all N particles being in the 
condensate, and none in the other M — 1 modes. Unfortunately, the Wigner function for 
an N particle state involves the Laguerre polynomial: 

2(—) N 

W N (a, a*) = exp(-2|a| 2 )L JV (4|a| 2 ), (83) 

and for N ^ this is a highly oscillatory function, as shown in Fig. [j], which has no 
probability interpretation. However, the cumulative integral 

Wff m (a,a*)= J d 2 a'6(\a'\-\a\)W N (a,a*) (84) 

behaves very like a step-function, and this means that the mean values of any smooth 
functions of a, such as the first few powers, will be very well approximated by using a 
Gaussian approximation of the form 

a^y^e 10 (85) 

where 9 has a uniform distribution on (0, 2n) and n is a real Gaussian random variable 
with 

(n) = TV +1/2 (86) 

var[n] = 1/4. (87) 

It should be borne in mind that in practice we are quite unable to even contemplate 
measuring any moments of a(x) higher than the fourth, and this form does correctly 
give all moments up to and including the fourth order. 

ii) The remaining modes would contain no particles, and thus for these the Wigner function 
( |83| ) is to be evaluated with N = 0, which is then is positive and Gaussian. 

iii) This means that we set 

a(x, 0) = v/ne^x, 0) + ^ £ fe (x, 0)a k (88) 

k 

= v^e i9 £(x,0) + y(x,0). (89) 
Here a k are complex Gaussian random variables, independent of each other, such that 

(o fc ) = (a* k ) = (90) 

(at) = (af) = (91) 

(a* k a k ) = 1/2. (92) 
It follows that that (p(x, 0) is a Gaussian random function with statistics given by 

(p(x,0)) =0 (93) 

(^(x,0My,0)) =0 (94) 

(^*(x,0) V *(y,0))=0 (95) 

(¥>*(x,0My,0)) =ip fl (x,y)-r(x,0)£(y,0)] (96) 

iv) The statistics of the initial random condensate wavefunction are Gaussian, with the non- 
zero means and correlations given by 

(a(x,0)) =£(x,0) (97) 

(a*(x,0)a(y,0))=iVr(x,0)e(y,0) + i<5 iJ (x,y) (98) 
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v) If there are M wavefunctions y n included in the summations (^TJ), and hence M— 1 terms 
in the summation over k in (p8|), then the mean number of particles Nw corresponding 
to the ensemble in (p7p8h is given by 



r M 
N w = I d 3 x(|a(x,0)| 2 ) =N+— (99) 

This can be seen correspond to N real particles, plus the effect of "half a quantum" per 
mode for the AI modes included. 

vi) If one runs a simulation for a time, and wishes to extract condensate and non-condensate 
number from the resulting ensemble, then the condensate wavefunction is 

*(x,t) = (a(x,i)) (100) 

and the condensate number is 

N c = J d 3 x\V(x,t)\ 2 . (101) 
The number of particles not in the condensate is 

/M 
d 3 x(\a(x,t)\ 2 )-N c -— : (102) 

that is, the "vacuum" particles must be subtracted. 

It is clear that the "vacuum particles" will contribute to the evolution by means of the nonlinear 
mixing arising from the particle interactions, and that their effect is quite possibly cutoff 
dependent. This cutoff dependence must disappear when one includes the full coupling to the 
noncondensate band by restoring the damping and noise terms in the second and third lines of 

3.4.2. Treatment of a non-pure condensate It is not difficult in principle to include a 
condensate and its quasiparticles, and this was done in the original treatment One simply 
expresses the operators afe in terms of the quasiparticle operators (3 k for the system in the 
Bogoliubov theory, and then then one assigns to these the appropriate mean values instead of 

(Pk) = (fit) = ( 103 ) 

Wi) - (K 2 ) = o (104) 

((3l(3 k ) =n k + 1/2. (105) 

Of course there will be considerable effort involved in solving the appropriate Bogoliubov- 
de Gennes equations in order to even construct the operators and their corresponding 
wavefunctions. However, this can be made much easier by using the method of Sinatra et 
al as explained in their later paper [[37]]. 



3.4.3. The choice of Er It is tempting to choose Er to be so large that the non-condensate- 
band phase space density F(H, x) is negligible, and all of the noise and damping terms in ( p3| ) 
can be neglected. However, this may not be possible, since this would mean that the higher 
energy parts of the non-condensate band would also need to have negligible occupation, and 
in this case the neglect of third order derivative terms would certainly not be permissible. Thus 
it turns out that the principal criterion for the choice of Er should be that there is significant 
occupation of all quantum states up to this level. 
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4. Approximate phenomenological equations 

The methods outlined above are quite complicated, and do not entirely answer the need for 
a simple set of equations which can give a quick estimate of the effects of being studied. 
Therefore we will consider in here a very simplified equation obtained by making some drastic 
approximations to the already approximate methodology we have developed. 



4.1. Mean value equations 



Let us therefore neglect the terms involving £ ,± , for the reasons noted above in Sect. 2.3.2 and 
consider the resulting equations of motion for the mean values: 

0(x,t) =T*{0(x,i)po(t)}. (106) 
nc(x,t) = Tr {0t( Xj t)0( Xj t) pc (t)} , (107) 

j c (x,t) =Tr(^ [0t(x,i)V0(x,t)- V[^( x ,t)]0(x,i)] p c (t) 

(108) 

The resulting equation includes a modified Gross-Pitaevskii equation, and a local growth 
equation. We include the relationship between the backward and forward rates G , which 
arises from the definitions ( j31p2| ) 

G<->[x,e c (x,t)] = e ^GW[x,ec(x,t)] (109) 
= J ||^V 2 0(x,i) - 7 T (x)0(x,t) - 2 W > NC (x,^(x,t) + U (0t( x )0(x) 2 ) 
+ {G( +) [x,e c (x,i)] (l-e e ° ( "T" M ) - M[x]} 0(x, t), (110) 
= V • J(x,t) + 2G«[x,e c (x ) i)] {(l - e S ^= Jt )n a ^t) + l} (111) 

Notice that 

i The terms involving G ± give local growth or decay of (f> and he, and they wall also cause 
some dephasing, as was analysed in QKIII and QKIV. When the forward and backward 
rates balance. 



ii The spontaneous term — the +1 in (1 1 1 ) — can initiate the condensate growth, while the 
term proportional to he gives the difference between stimulated growth and decay, which 
occur only with nonzero fie ■ 

iii The coefficient M[x] has no effect on the density he — it is a pure dephasing term of 



the kind well known in quantum optics and, as can bee seen from (109), it makes the 
amplitude of the coherent component <j> decay. Since it has no effect on the total number 
he, this means that the the effect is to transfer particles from the coherent component 
into the incoherent or thermalized component of the condensate band. 

iv) This kind of equation cannot be expected to give a good description of the full condensate 
growth process, from no atoms in the condensate up to a fully developed condensate 
in equilibrium with a thermal vapour. As shown in [[HJ ^4| ^9|| the full growth 
theory requires a much more detailed description of the kinetics of the thermal cloud, 
and one cannot simply assume the thermal cloud is always in equilibrium at a definite 



temperature and chemical potential. However, in [ 2 1 1 we noted that such a simple picture 



is reasonably valid once the condensate is reasonably large, say 50% of its final value, 
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and in [^5[ it was noted that the growth process is indeed well described as being at 
definite temperature, but that the chemical potential one should use is an "effective 
chemical potential", evaluated from the lower lying energy levels of the thermal vapour. 
Thus our description is probably a reasonable description of the process of matter 
exchange between condensate and thermal cloud in a situation not unreasonably far from 
equilibrium. 



4.2. The phenomenological growth equation 

The hydrodynamic approximation in practice should be defined in terms of the mean 
wavefunction 0(x, t), so that, to the extent that we can neglect derivatives of nc(x, t), we 
can write 

ec (x,t)=-»^M (112) 
at 

~ iH 9K ^ t] (113) 



0(x, t) dt 
We further make three further approximations: 



i) The argument of the exponent in (109) is sufficiently small to use e x m 1 + x. This is in 
practice almost always true. 

ii) We factorize all averages of products of operators. 

iii) We neglect completely the term M[x]<^(x, t) in (110), since it is expected to be small, 
and its effect is merely to cause a small change in fi. 



Doing these, we can replace (1 10) by 

<9<^(x, t) i ( ft 2 2 



dt h 1 2m 



V 2 0(x, t) - V r r(x)0(x, t) - 2wp NC (x, t)0(x, t)-u\ 



^<-^(x,i)-i^%^H. 



m 

where, for consistency with the notation of our previous papers, we have written 

W+ = G<+>[x,eo(x > i)]. (115) 

The function G + is slowly varying in space and time, so that W + can probably be 
approximated by a constant, since the most essential time- and space-dependence has been 
included in the factor fj,<fi(x,t) — ih,d(f>(x,t)/dt. The correct equilibrium results from 
the vanishing of the coefficent of W + , which ensures that the time dependence of the 
wavefunction is cxp(—ifj,t/h) — the remainng terms in the equation then reduce to the time 
independent Gross-Pitaevskii equation with chemical potential fi. Thus we get a stationary 
condensate wavefunction with the same chemical potential as the noncondensate. 



5. Application to vortex array stabilisation 

It is numerically known [ ^0[ [l2[ ^] that merely stirring a condensate described by 
the Gross-Pitaevskii equation produces vortices, but these do not stabilise into an regular 
array of vortices — as is experimentally observed [fS| [|5[ [ifj — merely as the result of the 
ongoing progress of the solution of the time-dependent Gross-Pitaevskii equation. To create 
a vortex array numerically, one solves the time-dependent Gross-Pitaevskii equation with 
an imaginary time and an appropriately time-dependent added chemical potential, while 
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constantly renormalizing the wavefunction so as to maintain a constant number of particles. 
This process finds a state with a energy minimum, but of course the method is purely an 
artefact, and does not represent the true underlying physics. 



However the form of the phenomenological growth equation (114) already includes what 
amounts to an imaginary time term, as well as a real time term. We first add an angular 
momentum term to transform to the rotating frame, and then use a simple form for W + as 
given in [ |l8| ] 

, Am(akT) 2 

W+(N)^g K —^-L. (116) 

nn 

where a is the scattering length for the nonlinear interaction term (so that u = Airah 2 /m), 
k is Boltzmann's constant, and T the temperature of the noncondensate). We have included 
the correction factor g w 3 to give an approximate match with the more detailed treatment, 
as suggested in [pip. The physics of this situation is the growth of a condensate in a frame 
rotating with angular velocity fl about the z-axis from a vapour cloud which is itself stationary 
in the rotating frame; thus in the laboratory frame, this is vortex nucleation from a rotating 
vapour cloud, such as has been experimentally implemented by the JILA group ffify with no 
rotating trap potential. 

Thus we find that we get 

dtp h 2 2 , . , ,2 

Here 



(*-7)&iir= -^-VV + ^r(x)^ + u|vr^-fii»^ + *7^NC^. (H7) 
at 2m 



4mga 2 kT 

1^ \ 2 ■ (H8) 

and using the values a m 10 _8 m and g = 3, we find the factor 7 r* 0.01. 

As noted in the previous section, the stationary solution of this eqaution will come from 
the equality of the coefficients of 7 on both sides, ensuring that the time dependence of the 
wavefunction in the rotating frame is exp(— i[/,nct/fi)> leaving the wavefunction to satisfy the 
stationary Gross-Pitaevskii equation (modified by the angular momentum term) with chemical 
potential /l«nc- This stationary solution will be a vortex lattice when it exists. 

More generally, we can consider a trap rotating with angular velocity £1 nucleating from 
a cloud rotating with angular velocity a. In this case we get, in the frame rotating with the 
trap, 

dtb h 2 2 

( z - 7) h ^7 = - tt V ^ + Vt(*)iI> + uM V - fiL,V + Z T + (o - fl)L z } ip. 
at 2m 1 1 

(119) 



The first equation < \l 17[ ) is similar to that presented in the paper of Tsubota et al. [[48[|, but 
there are significant differences: 

i) There is no physical justification given in [Q, and therefore their results, though 
qualitatively attractive, cannot be accepted as an explanation of the vortex lattice 
stabilization process. In contrast, our reasoning shows their choice of 7 = 0.03 is of 
the same order of magnitude as that expected from quantum kinetic theory. 

ii) They have a real term — fj/ip on the RHS, whereas this equation has an imaginary term 

iii) Their /i is adjusted with time to preserve the number of atoms in the condensate. Our 
/Unc is the physical chemical potential of the surrounding vapour, for which we have 
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no particular model at the moment. Even in the experiments [ |45| , |46| ] where an almost 
pure condensate is stirred by a rotating trap potential, it is known that as the stirring 
starts, there is considerable heating, and then at the end possibly 20% of the atoms form 
a thermalized rotating cloud, which appears to rotate at a lower speed than the trap. 



Thus we can expect that (119) should give a semiquantitative description of the actual 
process of vortex formation and decay when appropriate (possibly time-dependent) 
values of /knc and a are chosen. (However, we note that recent experimental work 
[ p^ l suggests that the stabilization process is only weakly dependent on temperature, 
which is surprising, since the dissipation in our equation is represented by 7, which is 
proportional to temperature. A resolution of this anomaly must await our more detailed 
calculations which are presently in progress.) 

The basic mechanism arises from the spatially dependent local chemical potential of 
a random lattice of vortices. The equilibrium lattice is characterized by a uniform 
local chemical potential, which must balance that of any surrounding vapour. The 
phenomenological growth equation includes this irreversible process. 

We can expect a model with fixed /inc to explain the actual observed dissipation of the 
vibration of vortices which is observed, but the full process of formation and heating could 
well be a formidable task, involving both kinetic and GP methodologies. Preliminary 
calculations have shown that qualitatively, the results are very similar to those of Tsubota et 
al. — that is, the formation and stabilization of a vortex lattice are observed in very much 
the same way. The numerical techniques necessary for both equations are almost identical 

6. Application to hydrogen condensate system 

6.1. Nonlinear losses 



The experiments on hydrogen 1 50, 51], |52j give rise to a system with a very large proportion of 
non-condensed hydrogen, which feeds a relatively small condensate as the condensate atoms 
themselves leave the condensate because of two-body dipolar relaxation. The data indicate 
that the condensate appears to have a density profile which matches that of a solution of 
the Gross-Pitaevskii equation. In contrast, the non-condensed vapour has a profile which 
does not match that of a zero chemical potential Bose-Einstein distribution, but instead 
seems to have significantly more population at lower energies than was expected from the 
zero chemical potential Bose-Einstein distribution. From the point of view of the quantum 
kinetic theory of condensate growth processes the zero chemical potential model chosen is 
not very appropriate, since any condensate has a positive chemical potential, and that of the 
noncondensed atoms must be even higher in order to maintain a net inflow of atoms into the 
condensate to replace those which are continually leaving by two-body processes. Thus, to 
the extent that a condensate in this situation can indeed be described by a chemical potential 
fic> one expects that the main body of the noncondensate would have a chemical potential 
Mnc > Mc> an d that there would be a transition region for the lower energy particles in the 
noncondensate. More precisely, what one expects is a distribution function f(E) of the form 

f{E) = cM{E-lAE))/kT]-l (120) 

where the energy-dependent chemical potential n(E) takes the value /i^c for higher E, but 
approaches hq as E — ► — the situation is illustrated in figure ||. 
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Figure 2. Illustration of the interpolation between a phase-space distribution form at high 
chemical potential ^nc at high energies, and the phase-space distribution at low energies, 
whose chemical potential /icis the same as that of the condensate 



6.1.1. Phenomenological growth equation including two-body losses We are most interested 
in the case of dynamic equilibrium, where the stationary state of the condensate is maintained 
by the flow of the atoms from the high chemical potential /znc through the vapour at chemical 
potential /Ltc and thence through the process of two body dipolar relaxation to what amounts to 
chemical potential of — oo, since this process is irreversible. In this case we have to consider 
three principal processes, namely 

i) The linear growth process, by which an atom leaves the vapour and enters the 
condensate. 

ii) The linear loss process, by which an atom returns to the vapour from the condensate. 

iii) The nonlinear loss process, by which pairs of atoms leave the condensate by dipolar 
relaxation. 

The first two of these terms are already included in the phenomenological growth 
equation (114), The nonlinear loss is obviously modelled by a term proportional to — \if;\ 2 ^j, 



leading to a phenomenological growth-loss equation for the hydrogen system 



h 2 _, . , |2 



ihijj = \7 2 ip + Vt(x)iI> + u\tp\ ip + ihiW 

2m ' ' 



kT 



R\il>\ 2 ^}. (121) 



It is a major advantage of the phenomenological approach which we have developed that 
the nonlinear losses can be accommodated so straightforwardly, with an effectively complex 
scattering length. The simplicity of the modified Gross-Pitaevskii equation allows quantitative 
investigation of questions which might well be difficult even to formulate, if one had to 
proceed directly from first principles without the phenomenological theory as an intermediate 
step. Yet among such questions are those with potentially dramatic physical consequences, 
such as the one we now pursue. 



6.2. Instability of the Thomas-Fermi stationary state under nonlinear loss 

Under currently typical experimental conditions, the timescales of condensate growth and 
decay, including two-body loss, are much longer than those related to the harmonic trapping 
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potential, the mean field energy fj,c, or the non-condensate chemical potential ^nc- To take 
advantage of this small ratio of frequency scales, it is usual to introduce length and time scales 
related to fic (the healing length at peak density, and the mean field time scale); but since in 
the hydrogen context we want to see how the condensate maintains dynamic equilibrium in the 
presence of a much larger thermal cloud, we will consider //no as the fixed parameter, which 
determines /ic- Since in dynamic equilibrium the two chemical potentials are comparable, 
this difference is not important. So we will convert to dimensionless variables 

t = (122) 



x = - — x (123) 

n 

V(Z) = Mnc^W ( 124 > 

u = /%c u • (125) 
This also motivates the dimensionless parameters 

c = 2W + ^ (126) 
uW+ 



so that we can re-write ( 121 ) as 

1 ~ ~ IP 

iip = - - + u\ip\ 2 i> + V{Z)il> + —(ip -tip- uA\ip\ 2 i>) . (128) 

Since W + is on the order of the Boltzmann scattering rate ap tc VT for the thermal cloud, for 
the hydrogen condensate at temperature approximately 50 /iK this gives e ~ 10~ 6 . And with 
the published values for the two-body loss rates in this system, the dimensionless parameter 
A is on the order of 10. The analysis below suggests that this is quite large enough to raise 
serious questions about the stability of the Thomas-Fermi stationary state; and it is possible 
that the large occupations in the lower quasiparticle levels may mean that A ~ 10 is an 
underestimate. 

In the limit where e — > 0, we obviously obtain the standard Gross-Pitaevskii equation. 
Since the trapping potential varies very slowly on the healing length scale, in this limit the 
Thomas-Fermi approximation is excellent. If e is of order unity or larger, the Thomas-Fermi 
stationary solution breaks down badly, and it may not even be possible to find a stationary 
solution to (128). But in the experimental regime of small e, there are indeed stationary 
solutions in which the density profile is very close to the standard Thomas-Fermi form. 
However, unless A is small enough, they are dynamically unstable. By itself, the two-body 
loss term drives the total number of particles towards an equilibrium value. And it even tends 
to correct against collective excitations of the steady state condensate. In bulk, it is easy to 
show that this corrective effect strongly stabilizes against perturbations; but as we show below, 
in a harmonic trap it overcorrects, and the collective excitations grow in amplitude. Numerical 
integration reveals that these excitations continue in 'boom-and-bust' cycles (where 'bust' 
means a collapse of central density to a fraction of its maximum value); but in this regime, the 
assumptions that justify the phenomenological mean field theory may well be breaking down. 
In an actual experiment, it is not clear whether we should expect coherent but non-stationary 
states like those assumed in our theory, or some kind of quasi-condensate with degraded phase 
coherence, or merely accelerated decay of the condensate. 
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6.2.7. Stationary states We begin by extending the Thomas-Fermi stationary state to first 
order in e. Since it is still true that V is a slowly varying function, we are in the hydrodynamic 
limit, and we can approximate (128) by the hydrodynamic equations for ip = y/pe te and 
v = V6*, namely 

p = - V ■ (pVO) + e(l + - ukp)p (129) 

0= -\\V6\ 2 -up-V(Z) , (130) 

where we drop a term — ep/(4p) in the second equation. (We drop this term because, even 
when below we allow time-dependent density perturbations around the stationary state, this 
term will be of order e 2 . It does become large as we approach the Thomas-Fermi surface, but 
this only leads to a small correction to the usual boundary layer theory that must be applied 
to match the hydrodynamic approximation to the near-surface region.) A stationary solution 
allows = —pc/PNC as the only time dependence in ijj, giving 

V • (pv) = e(1 - nc/nxa ~ uAp)p (131) 
^ ^+(^V?+(^Vf + 0(e 2 ) (132) 

for cylindrical co-ordinates and an axisymmetric harmonic trap. So to first order in e, we 
have the familiar Thomas-Fermi density profile for an axially symmetric harmonic trap: 



1 

P = - 

u 



r, 2 ; 2 

U) , Z 



(133) 



where we introduce the dimensionless trap frequencies cDj= (Jj/pwc an d condensate 
chemical potential p= p/p-^c- But there is now also a velocity field W9 of order e. 

Since a very small velocity field is difficult to observe, we are mainly interested in the 
density profile. Actually finding v to first order is therefore only indirectly necessary (it fixes 
pc)\ but it is not so hard. Assume an ansatz of the form 
d9 

— = er(a r + b rr r 2 + b rz z 2 ) (134) 
or 

d9 

— =ez{a z + b zr f 2 + b zz z 2 ) (135) 
oz 

for unknown constants aj,bjk- Since W must be irrotational, b rz = b zr . This leads to 
six equations for only five free parameters in the ansatz, and so forces a specific value on the 
heretofore undetermined term in p, namely the condensate chemical potential p (which makes 
sense, since we do expect the growth and loss terms to determine pc). The result for the 
axisymmetric trap is 

80__ e (7Cj 2 + 4u 2 )(p - cD 2 f 2 / 2) - f Co 2 z cu 2 P = 

df 7 5a; 2 + 6ui, 

80_ e (4£ 2 + 7u 2 z ){~p - Cj 2 J 2 /2 ) - ^Cup 2 f 2 

dz 7 5(D 2 + 6ij 2 

= IT|a ' (138) 

The uniqueness of this ansatz solution can be proved straightforwardly in the ID limit, and 
also for spherical symmetry; we conjecture that it is more generally unique. 

The results of (136-138) are somewhat complicated, but they make excellent physical 
sense. Our expectation that the nonlinear losses will lower pc below p^c is borne out. And 



f (136) 
-2 (137) 
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Figure 3. Hydrodynamic flow lines in the presence of weak two-body losses, for 
harmonically trapped axisymmetric condensate with aspect ratio 2. The flow lin es are a 
seen in a section parallel to the symmetry axis. Because of the symmetry of Eqn. (|136Ul37 
under r <-> z, this Figure is identical whether the condensate is oblate or prolate. 



the velocity field is inward from the edges of the condensate towards the core. Since the loss 
rate scales faster with the density than the growth rate, the periphery of the condensate receives 
from the thermal cloud more particles than are needed to maintain the local condensate 
density, and so it is able to donate a flux to the central regions, where loss rates exceed growth. 
The velocity field becomes radial and proportional to p in the limit of spherical symmetry, but 
otherwise it has no simple description; see Fig|| for an illustrative example. In the limit 
uj z <C ui r , we find d z becoming independent of r, and the same statement with r and z 
interchanged is also true, so we do obtain ID and 2D limits where we expect them. 

And it is always true that Vf? • Vp = at the edge of the Thomas-Fermi cloud, which 
is required in the hydrodynamic approximation. (Density vanishing outside the TF surface 
means no flux can cross, since gain and loss terms both vanish at zero density). 

So with the addition of the slow flux from periphery to core, the Thomas-Fermi solution 
is essentially maintained in the presence of slow two-body losses. But is the Thomas-Fermi 
solution stable? 



6.2.2. Collective excitations To answer this question we need to find the frequencies of 
linear collective excitations, by linearizing (129,130) about the stationary state (131 132). 
Using a method which can be dignified as multiple scale analysis, but which is essentially the 
same as time-independent perturbation theory in ordinary quantum mechanics, we compute 
the frequencies of the collective modes to first order in e. The first order corrections are in 
general imaginary (as is not surprising given that our perturbations are imaginary). What may 
be surprising is that if A exceeds a threshold of order unity, the amplitudes of some collective 
excitations grow, indicating that the two-body losses cause an over-corrective instability. 
The linearized hydrodynamic equations, to first order in e, are 

5p= - V ■ (pVSO + Sp V0) +e(l - p, - 2Aup)5p - pSO 

SO = -u8p-V0-V50 . (139) 

We can combine these into a single second order equation for 80 (again dropping higher order 
terms in e): 

80 = uV • (pVSO) - 80 V 2 - 2V6> ■ VSO - e[l 



The idea now is to look for 



80{±), 

eSOi + .... At zeroth order we have 
V- l/2-fr 2 - 



ng<w 



- /i - (2A - 1)^]^ . 
and expand Q = Qq 



(140) 
., and 

(141) 
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Then using the orthonormality of the eigenfunctions of the zeroth order RHS, we obtain as in 
the time-independent perturbation theory of quantum mechanics 

-2O Oi / d 3 f-56l= -ifi / d 3 f[l - p - (2A - l)up]8dl 
Jtf Jtf 

10 



where the integral is over the Thomas-Fermi cloud. (The terms on the RHS with V0 in them 
add up to V • (69q V#), which by Stokes' Theorem is a surface integral that, as we mentioned 
above, vanishes, since the normal to the surface of the TF cloud is oc Vp.) 



If A < 0.7, then our stationary solution is stable because the integrand in ( 142 ) is positive 
throughout the volume of integration; but if A > 0.7, then the integrand is negative in a 
volume around origin, and so the integral might perhaps be negative at least for some modes. 
To investigate further we must learn the 50$ for the various collective modes. 



Solving (141) for the hydrodynamic modes of a harmonic trap can be reduced to 
diagonalizing finite matrices using a Frobenius series approach. In an axially symmetric trap 
it is not hard to find modes up to third order polynomials analytically and exactly, and with 
spherical symmetry we can find all the modes. In an extremely prolate (cigar-shaped) trap, 
we can also find all the modes to zeroth order in (tu z /ui r ) 2 . We will proceed here with the 
spherical calculation. The extremely prolate case is analysed in the Appendix, and here we 
will only quote its results. 



6.2.3. Spherical trap In the spherical case it is convenient to solve (141 ) in spherical polar 
co-ordinates, with the angular dependence of 89q obviously being a spherical harmonic Yi m . 



With oj. 



ui, define the rescaled spherical radius 



and write 59q = Yi m (d, 4>)f n l( r ) s ° that (141 ) becomes 



2ng f 

~ 9 Jnl 



[l(l + l)(l-r 2 ) + r- 2 d r (r 4 -r 2 )d r ]f nl . 



(143) 



(144) 



(The reason for the additional subscript n on f n i will appear shortly.) Taking a Frobenius 
ansatz /„; = akni r fc+2 for the radial dependence then yields the recursion relation 

[(k + s + 2) (k + s + 3)-l(l + l)\ 

= [{k + s)(k + s + 3) - 1(1 + 1) - 2n 2 /Co 2 ] a kn i. (145) 

Regularity at the origin then forces s = I, and further requires that akni vanish for all k with 
parity opposite to I. Convergence as we approach the Thomas-Fermi surface r = 1 forces the 
sequence of non-zero akni to terminate at some k = n — I > 0. This introduces the quantum 
number n, which must have the same parity as I, and be greater than or equal to it. We recover 
the familiar result 

2^ = Co 2 [n(n + 3) - 1(1 + 1)] ; (146) 

the point of repeating the calculation to this point has been to recall the recursion relation 
(145), from which we can evaluate the RHS of (142). 



In our rescaled spherical co-ordinates the equation for Oj becomes 



sphere 



ip 

~2 



10 

y 



A - (1 -2A) 



(147) 
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We can evaluate the ratio of moments in this expression by using the orthogonality under the 
weight r 2 of the eigenfunctions of ( 144) of different n, together with a fact we can extract from 
the recursion relation (145). Since the /„; are polynomials of rank n and definite parity, it 
follows that 



fn+2,l = J]] 0>k,n+2,l 
k=0 



k+l 



dn-l,nl 



r fnl 



(I 



for some constants Cj 
r 2 fn 



n+2-Ln+2,l 

Q>n — l,nl 
n-2 

This then means that 



n-l,n+2,l 
0~n+2-l,n+2d 



a n-2-Lnl 



&n-2~l.nl 



a 



n—l,nl 



n-l,n+2,l 
0-n+2-l,n+2,l 



fnl + > , Cjfji, 



so that by orthogonality we have 



Jo ^ r 1-4 fnl _ a n-2-l,nl 



i 



0"n—l,n+2,l 
0-n+2-l,n+2d 



( n+ l)(„+|) 



where the evaluation follows straightforwardly from (145). 
This immediately yields our desired result 

1-®A+(2A-1)- 



sphere 



Ifi 

T 



(n+|)(n+|) 



(148) 



(149) 



(150) 



(151) 



(152) 



(153) 



(154) 



(155) 



^From this we can see that the so-called 'surface' modes, with I = n, are always stable. (The 
case n = I = is not an excitation, but simply the stationary state, so this case does not 
count.) Since if A < 1/2 then Imf^ phcrc will be positive because 1 — 6A/7>1 — 2 A, 
we need only look for instabilities in cases where 2 A — 1 > 0. In these cases it is clear that 
instability can only occur for A > 7/6, in which case it occurs first for smaller / and larger 
n. Equation (155) indicates that for A > 7/6 there will always be instabilities at sufficiently 
high n, but the hydrodynamic approximation on which the equation is based will break down 
for n larger than some limit that depends on /xc and u>, and so in general the actual instability 
threshold will be somewhat higher. But for I = 0, all modes with n > 2 will be unstable 
if A exceeds only 77/64, and so the instability threshold is not actually very sensitive to 
post-hydrodynamic corrections. 



6.2.4. Extremely prolate axisymmetric trap Instability is not a pathology of an exactly 
spherical harmonic trap: it occurs also in a quasi- ID trap (for which all of our calculations 
can be repeated quite simply), and in an extremely prolate but hydrodynamic ally three- 
dimensional trap (which is not the same thing). For the extremely prolate axisymmetric 
trap the result obtained in the Appendix is 



pro 



Ifi 

~2 



iA+(l-2A) 



p 2 + p{2n + 3) + (re + 2)(2n + 1) 
(2p + 2n+ l)(2p + 2n + 5) 



1 



re(re + 2) 



(156) 
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where |m|, n,p are whole numbers. The azimuthal quantum number m must have the same 
parity as n, and \m\ < n. In this case, all of the axial modes (n = m = 0) are stable even 
if A — > oo. And the dipole modes (n = m — 0,p — 1 and |m| = n = = 0) are not 
only stable, but have Ox independent of A. (This is also true in the spherical case, and indeed 
for all harmonic traps, because the dipole modes merely translate the entire condensate, and 
so do not disturb the local balance between growth, loss, and flux.) Raising |m| and p 
tends to stabilize; the most unstable modes will have p = m = 0. For these modes we can 
recognize that a dynamical instability occurs for (6n + 8)A > (7n + 14), which will occur 
at high ?i for A > 7/6, just as in the spherical case. The same caution applies, that for n 
too large the hydrodynamic approximation breaks down. But for A > 7/5, all modes with 
p = m = 0, n > 2 will be unstable, and so again post-hydrodynamic effects cannot shift the 
instability threshold very much. In the hydrogen condensate experiment at MIT, where the 
trap has an aspect ratio of 400, the extremely prolate limit definitely obtains, and with A ~ 10 
there are clearly a lot of unstable modes. We conclude that instability of the Thomas-Fermi 
stationary state is a general phenomenon if there are sufficiently strong two-body losses. The 
important point here is that 'sufficiently strong' means only that the two-body rate constant 
be slightly greater than the growth rate constant. 

6.3. Implications 

Since the MIT hydrogen condensate experiments are well above the instability threshold, 
our phenomenological growth equation seems to indicate that the quasi-steady state of these 
condensates cannot be the Thomas-Fermi density profile with perfect phase coherence, which 
all other condensates exhibit well. The problem is not so much that the Thomas-Fermi state 
should be oscillating unstably, but that the system cannot be expected to settle down to an 
unstable state in the first place. 

How seriously should we take this conclusion, which is after all derived from 
approximate solutions to a phenomenological theory? We would argue that, even if the theory 
might not be quantitatively precise, the overcorrective instability which it predicts under two- 
body losses is a simple physical phenomenon. When loss and gain processes scale differently 
with density, and density is inhomogeneous, particle flow is needed to maintain a steady state. 
In effect the condensate velocity field is a control system that tries to maintain a constant 
density profile. But moving atoms have momentum, and so the problem of overcorrection 
can obviously arise in this control system; whether this leads to instability is a detailed 
matter of length and time scales. It is possible that our approximate theory exaggerates the 
overcorrection problem, or that it may be counteracted by factors not included in the theory; 
but it is a genuinely physical possibility, and not a mere artifact of our phenomenological 
approach. 

Among factors neglected in this Section, which might tend to suppress instabilities in 
real systems, is the thermal component of the velocity field. In the hydrogen experiments the 
thermal fluctuations in the condensate velocity field can be estimated to be much larger than 
the systematic flow which compensates for the two-body losses, and so these should really 
be taken into account, by solving the full, stochastic version of our modified Gross-Pitaevskii 
equation. 

It is a great merit of the phenomenological growth equation that it allows explicit 
calculation of the instability question, which would be difficult even to formulate directly 
from first principles. But given the experimental data showing a high-density component of 
ultracold hydrogen, whose density profile seems consistent with those of condensates in the 
alkali vapours, the phenomenological growth equation cannot be correct when it concludes 
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that the Thomas-Fermi ground state is impossible in hydrogen. The results of this Section 
do indicate, however, that the balance between thermal fluctuations and dissipation (loss) 
must be much less trivial in hydrogen than in the other ultracold Bose gases. Thus in this 
Section we have provided evidence that further development of non-equilibrium theory for 
cold, dilute, trapped bosons will provide insights into qualitatively new regimes, and not 
just small corrections to familiar results. At the same time the stochastic Gross-Pitaevskii 
approach has proven itself as a workable tool which can, without requiring unreasonable 
effort, yield definite answers to nontrivial physical questions. 



7. Conclusion 



In this paper we have given indications of how a stochastic Gross-Pitaevskii equation should 
be defined and used in practice. There are two major issues to be resolved in this formulation 

a) What is the appropriate definition of the "condensate wavefunction" i/j{x., t) for which 
the stochastic Gross-Pitaevskii equation should provide an equation of motion? Our 
conclusion is that a multimode Wigner function definition is the most appropriate, 
and although this cannot be done exactly, or even for all possible quantum states, the 
approximations and restrictions which are thus implicit in this definition appear to us to 
be unlikely to cause in serious problems when applied to practical situations. The most 
irksome complication is the necessity to include "vacuum noise" in the initial conditions, 
basically in ord er to make sure the Heisenberg uncertainty principle is not violated, as 
detailed in Sect. 33. However, at all but the very lowest temperatures this is unlikely to 
be an issue — the influence of invalid initial conditions will very soon be eliminated by 
the noise induced by the thermalized atoms. 

b) The formulation a quantum kinetic master equation using the local formulation of energy 
and momentum conservation introduced by Zaremba et al jll] ] creates a very much 
simpler master equation than that of our earlier formulations, and its transformation 
into the stochastic Gross-Pitaevskii equation is then relatively painless. There are, 
however, many technical details involving the best way of quantitatively specifying the 
noise coefficients which are still in need of refinement. We do not see these as being 
too important quantitatively in practical situations, but they will eventually need to be 
attended to more precisely than we have done here. 

The major result is that, subject to a number of caveats, including the correct choice of initial 
conditions, we can modify the Gross-Pitaevskii equation to include damping as a result of 
exchange with non-condensed atoms, in order to obtain a semiphenomenological description 
of the interacting systems which includes all of the quantum effects in at least an approximate 
form. 

The development of the phenomenological growth equation, is similar to one presented 
by Williams and Griffin [^6||, is the most useful immediate application of this work. This 
equation can be seen as a very simplified description of the interaction of a condensate and a 
thermal cloud, which includes the major processes in a consistent but only semiquantitative 
way. The formulae for vortex lattice nucleation and stabilization given in Sect.|| are new, and 
in work which we shall publish elsewhere, we will show how a rich variety of behaviours 
arise out of this simple and elegant formulation. 

We have also exploited this equation investigate the description of growth and loss in 
a hydrogen condensate, and have shown that there are hitherto unremarked hydrodynamic 
instabilities in the Thomas-Fermi stationary state of a cold condensate suffering two-body 
losses. We shall publish more detailed numerical studies of this system elsewhere. 
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Appendix: two body instability in extremely prolate traps 



Zeroth order modes and overlap integrals The first step is to introduce the rescaled co- 
ordinates 



R 



and define 20q 
becomes 



•qui, 



vf(R, Z) = (d f 



(157) 

f(R, Z)e un ^, so that our zeroth-order equation 

2 (1 



R-i)(l- 
+ Vdz(l 



R 2 



R l 



Z 2 )d R f - m 2 R~ 
- Z 2 )d z f 



R 2 



Z 2 )f 
(158) 



and our frequency correction becomes 



4 -A 
7 



(1 - 2A) 



Io dZ Jo 



RdR (1 - R 2 - Z 2 )f 



Stdztf 1 - 2 RdRf 2 



■(159) 



With an aspect ratio of 400 in the hydrogen experiment, we can certainly stop at 
zeroth order in r). To do this it is convenient to introduce yet a further new co-ordinate: 
R = s/l — Z 2 X. (This necessitates some care with expressing partial derivatives and 
measures in the new variables, but the results can be checked by comparison with exact 
solutions to the equations in the original variables, and we have done this for the ten simplest 
modes.) In these final variables X, Z we have 

- uf(X, Z) = (d x + A- X )(l - X 2 )d x f + m 2 (l - X- 2 ) 



1 1 1. 

~~2 



-A + (1-2A) 



J^- 2dx]( l-Z 2 )(l-X 2 )(dz + ^- 2 
tidZ(l-Z 2 ftiXdX(l-X 2 )f 



tidz(i-z*)tixdxp 

G(Z)F(X) + 0(ri), v = vq + i)vi 



dx)f (160) 
(161) 



It is therefore clear that f(X, Z) = G{Z)F(X) + 0(rf), v = v a + i)v x + 0(rf) 
- v F = (d x + - X 2 )d x F + m 2 (l - X- 2 )F . 

Writing F(X) = ^ akX k , we get the recursion relation 

[(k + 2) 2 - m 2 ]a k+2 = [-v + k(k + 2) - m 2 ]a k . 

Since the series can only converge if it terminates, we conclude that v Q = n(n + 



, with 

(162) 

(163) 
2) — m 2 for 



some integer n. Since F must not blow up at X = 0, we must eliminate negative powers of 
X; this requires that a k = for all k < \m\. This implies that, for a non-vanishing solution, 
\m\ < n; and since the recursion relation goes in steps of two, this further implies that all 
non-vanishing solutions have m the same parity as n. 
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We can use Q162J) to show that solutions F nm with different n are orthogonal in the 
interval X G [0, 1] with the weight X. Consequently, we can evaluate moments such as 
/ dX X 3 F 2 m / J dX XF 2 m using the recursion relation (163), by observing that 



X 2 F nj 



/] a n -k- 



iX k 



-F n 



+2,m 



a-n+2.n+2 

for some q. This implies that 



*n,n— 2 
ii .n 



a 



n+2,n+2 



qFi 



3p2 
n.ra 



JodXXFl 
We can now go back to dj~6 



(n + 2) 2 - m 2 
4(n + 2) 



4n 



ra(ra + 2) 



(164) 



.(165) 



to first order in 77, and by integrating over X with XF nm , 
proje ct out the equation for G(Z). (We have to integrate by parts several times, and apply 
(165) as well.) In fact we find a simpler equation for g[Z) = G(Z)(1 — Z 2 )~ n l 2 (which is 
reassuring since after all X n = R n (l - Z 2 )~ n / 2 , and so F nm (X)(l - Z 2 ) n / 2 g(Z) will be 
regular at R = 1, Z — > ±1 if g is regular at Z = ±1). The result is 

{l-Z 2 )g" -2{n + 2)Zg' 



(166) 



where a fairly complicated bunch of terms involving m, n and v\ has been absorbed into the 
eigenvalue £. Similar analysis to that above shows ^ = p(p + 2n + 3) for non-negative integer 
p, determining v\ and hence as a function of the mode indices m,n,p. We have thus found 
the complete hydrodynamic spectrum of an extremely elongated harmonic trap. 

By a similar analysis as already performed for the F nm above, we can also compute the 
second moment of g(Z) (which are orthogonal over the interval Z G [—1, 1] with the weight 
(1 — Z 2 ) n+1 ). So we also finally have the correction to the frequency to first order in e: 

(p,n\(l - Z 2 )\p,n) z (n,m\(l - X 2 )\n,m) 



2 

ijl 
~ ~2 



4 -A- 
7 



yA+(l 



(1-2A)- 



IX 



2A) 



(p,n\p,n)z (n,m\n,m)x 
p 2 + p(2n + 3) + (n + 2)(2n + 1) 



(2p + 2n+ l)(2p + 2n 



1 - 



(167) 
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